Coherence and Josephson oscillations between two tunnel-coupled one-dimensional 
atomic quasicondensates at finite temperature 
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We revisit the theory of tunnel-coupled atomic quasicondensates in double-well elongated traps 
at finite temperatures. Using the functional integral approach, we calculate the relative phase 
correlation function beyond the harmonic limit of small fluctuations of the relative phase and its 
conjugate relative-density variable. We show that the thermal fluctuations of the relative phase 
between the two quasicondensates decrease the frequency of Josephson oscillations and even wash 
out these oscillations for small values of the tunnel coupling. 

PACS numbers: 03.75. Kk,03.75.Lm,67.85. Jk 



I. INTRODUCTION 

Systems of ultracold bosonic atoms in two parallel 
atomic waveguides mutually coupled via quantum tun- 
neling (so-called extended bosonic Josephson junctions) 
have been a subject of intensive theoretical QHZI and ex- 
perimental Q studies. The finite spatial extension of 
these systems provides much richer physics compared to 
the case of a point-like bosonic Josephson junction 
The novel features arise due to the enhanced role of noise 
and correlations in low-dimensional ultracold atomic sys- 
tems. 

Before discussing the effects of tunneling, we recall the 
basic proper ties of a bosonic system in an isolated waveg- 
uide [1CI - I12T ]. This system is effectively one-dimensional 
(ID), if the interaction energy per atom (we assume in- 
teratomic repulsion characterized by the effective ID cou- 
pling strength g > 0) and the temperature are well be- 
low the spacing between the discrete energy levels of the 
potential of tight radial confinement. In this case quan- 
tum degeneracy does not lead to establishment of the 
long-range order; instead, atoms form a quasicondensate, 
i.e. a system describable by a macroscopic wave func- 
tion with strong phase fluctuations. The characteristic 
length of the phase coherence in a quasicondensate at fi- 
nite temperature T is A^ = 2h? n\Y) / [rak^T) , where m 
is the atomic mass, and n^D is the mean linear density 
of atoms fioj (we assume an infinite system; thermody- 
namic limit implies constant nio = N/L while both the 
atom number iV and the quantization length L tend to 
infinity). The power-law decrease of the single-particle 
correlation function takes place only at T = 0. 

If two waveguides are tunnel-coupled, the system is 
described by the generalized Hamiltonian 
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waveguide (j = 1,2), /i = hgniD — TiJ is the chemi- 
cal potential and 2 J is the tunnel splitting (in frequency 
units), i.e., the frequency interval between the two lowest 
eigenstates of the radial trapping Hamiltonian (the an- 
tisymmetric and symmetric superpositions of the single- 
atom states localized in either j = 1 or j = 2 wells of 
the double- well Hamiltonian). In this case the situation 
changes qualitatively: the tunnel coupling mutually locks 
phase fluctuations in the two quasicondensates [l| . Phase 
locking (as we shall quantify later, in Sec. |H| means that 
the distribution of the relative phase between the two 
quasicondensates becomes peaked around zero, while the 
local phase of an individual (j = 1 or 2) quasicondensate 
remains fully random (the phase-density representation 
for quasicondensates will be discussed in Sec. [II]). In the 
spatial correlation of the local relative phase between two 
quasicondensates a new length parameter appears [H [f| , 



lj = y/n/{AmJ). 



(2) 



where ipj is the atomic annihilation operator for the jth 



The length /,/ sets the scale of restoration of the 
inter-waveguide coherence due to finite tunnel-coupling 
strength J. The tunnel-coupling strength is usually es- 
timated from the single-particle energy (kinetic and po- 
tential) and the overlap in the potential barrier region of 
the wave functions for a particle localized in the 1st and 
2nd waveguide. However, it is also possible to take into 
account atomic interactions, see Ref. @ and references 
therein. 

Experimentally, the interwell coherence can be ob- 
served by releasing the two quasicondensates from the 
trap and measuring locally the contrast and t he p hase of 
their interference pattern after time of flight [l3l fl4|. 

Up to now, only the theory based on linearization of 
the Hamiltonian ((1) has been developed Jl land applied 
to the analysis of the experimental data [a, If . Our work 
is aimed to develop a model of the steady-state thermal 
noise in tunnel-coupled quasicondensates beyond the har- 
monic approximation as well as to quantify the influence 
of the thermal noise to the macroscopic coherent dynam- 
ics of the system (Josephson oscillations). 

Our paper is organized as follows. In Sec. |H]we sum- 
marize the harmonic approach of Ref. [l[ . Section Mil is 
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divided in two Subsections. In Subsection IIII Al we calcu- 
late the static correlation properties of our system beyond 
the harmonic approximation using the transfer operator 
technique in the classical limit. The condition for negli- 
gibility of the quantum noise is also derived. Our way to 
model equilibrium state by numerical simulation of the 
system's relaxation to the equilibrium after a quench is 
explained in Subsection IIII Bl Section IIVI deals with the 
noise-affected Joscphson oscillations. We derive analiti- 
cally the frequency of Josephson oscillations modified by 
the thermal noise in our extended system. We support 
our analysis by numerical simulations and also observe 
Josephson oscillations washing out as a result of thermal 
noise for low enough tunnel coupling. Section IVl contains 
our final remarks and conclusions. Explanations of the 
ways to derive the main equations of Sec. IIIII and Sec. 
IIVI are presented in Appendices [5] and |B] respectively. 



II. HARMONIC APPROXIMATION 

Following the standard procedure [llj, we represent 
our atomic field operators through the phase Oj(z) and 
density Pj(z) operators, obeying the commutation rela- 
tion [Oj(z), Pj'(z')] = —iS(z — z')Sj f, as 



i>j (z) = exp[i0,-(.z)]^/ 'pj(z), 



1,2. 



(3) 



A discussion of the way to introduce the phase operator 
for quasicondensates by coarse graining a lattice model 
on length scales containing sufficiently many atoms can 
be found in Ref. 11 1. The density operator can be rep- 
resented as Pj(z) = niD + Spj(z). Since for quantum 
gases with repulsive atomic interactions density fluctua- 
tions are suppressed, we can always consider the corre- 
sponding operator Spj as a small correction. However, 
the same is not always true for the phase fluctuations. 

Whitlock and Bouchoule [l[ from the very beginning 
assumed the phase fluctuations to be small and thus lin- 
earized the Hamiltonian Eq. ([I) reducing it to H m H\[ n , 



lin 




gSp a + TiJnmO 



Here the symmetric (s) and antisymmetric (a) variables 
are introduced via canonical transformation 

6p s (z) = Sp 1 (z) + Sp 2 (z), e s {z) = [e 1 (z) + e 2 (z)}/2, 

Sp^z) = [Sp^z) - 6p 2 {z)]/2, 6 & {z) = e 1 (z)~e 2 {z). 

Diagonalization of the Hamiltonian (j4|) is based on 
the Fourier transform <5/5 s(a) (z) = L" 1 / 2 J2k^o ^s{a,),k^ %kz , 



Q s (a)(z) = L 1/2 Efc#o 4(a), fce lfcz . The frequencies 
w s ( a )(fc) of the symmetric and antisymmetric modes with 
the momentum hk are given by the dispersion relations 
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Correlations in two tunnel-coupled quasicondensates 
are experimentally accessible via the two-point correla- 
tion function g\{z—z') — n^^.^\(z)^\(z')^ 2 (z)^)\[_z') :). 
Since the system described by the Hamiltonian ([T]) is 
translationally invariant, g 2 depends only on the differ- 
ence of the two co-ordinates. The symbol (: O :) denotes 
the average of the normal ordered (with respect to the 
atomic operators tpj, ?/>j) form of the operator O. In what 
follows, we omit the normal ordering notation, thus ne- 
glecting the atomic shot noise. 

Since the density fluctuations for |fc| < £ = 

h/ ^JmgniD = h/{mc) being the healing length, are sup- 
pressed by the atomic repulsion fiol [Tl| , the main contri- 
bution to this correlation function is given by the phase 
fluctuations, g 2 (z — z 1 ) ss (exp[i# a (z') — i0 a (z)]). 

The experimentally accessible length scale cannot be 
shorter than the optical resolution length Az opt . On this 
scale the shot noise yields the quantum uncertainty of the 
relative phase, coarse grained over the distance Az opt , of 
the order of 1/ ^/2niDAz opt . For Az opt ^ 3 pm and 
"■id ^ 30 /an" 1 the shot-noise induced phase uncer- 
tainty does not exceed 0.075 rad. This relatively small 
value can be always kept in mind when comparing the- 
oretical predictions to measurement results. However, 
for the sake of simplicity, in what follows we assume 
(: exp[ie a (z') - i6 a (z)} :) « (exp[i6 a (z') - *&(*)]} and 
so on. 

Another point related to the use of the fully classical 
approximation is the substitution of the Bose-Einstein 
statistics of the elementary excitations by its classical 
limit, 



fc B r 



exp[hu a (k) I (k B T)} - 1 huj a (k) ' 



(7) 



One obtains strong deviations from Eq. (JT]) for hoj a {k) > 
k^T, which corresponds, under typical experimental con- 
ditions, to the range of wave lengths shorter than Az opt , 
i.e., not resolvable optically. 

These considerations justify our method based on gen- 
uinely classical statistics. 

In the harmonic approximations fluctuations arc Gaus- 
sian, hence, (exp[i6> a (z') - i6 a (z)}) = cxp{-±([0 a (z') - 
9 a (z)) 2 )}. Expressing 8 a through creation and annihila- 
tion operators of the elementary excitations and calculat- 
ing thermal populations of the elementary modes using 
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Eq. (O, Whitlock and Bouchoule obtained [l[ 



ate the partition function [Io| 



^a(^)]) = exp 



2h 
At 



(1 



,— \z — z'\/lj 



Z = J VS Ps J V9 S J V5 Pa J V9 a exp[-H/(k B T)}, 



(8) 

From this expression we can see that tunnel coupling 
locks the relative phase between two quasicondcnsatcs. 
This locking means that the relative-phase correlation 
function ((SJ docs not decrease to zero, but even at 
\z— z'\ — > oo has a finite value, corresponding to (9 2 (z)} = 
21j/Xt- On the contrary, the phase correlations in 
each of the waveguides are (exp[i6j(z') — iOj(z)]) = 
(exp{i[6 s (z') ± \9 a {z') - 9 s {z) t \0 a {z)]}), the upper and 
lower signs corresponding to j = 1 and j — 2, respec- 
tively. We can evaluate them using the statistical inde- 
pendence of noise in the symmetric and antisymmetric 
modes. The result 



(10) 



where 
H = 



Ti 2 nm ( d0 s \ g 2 



% niD / d6 
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gSp 2 + 2hJniY)(l — cos0 a ) 
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is the Hamiltonian ([T]) expressed through the classi- 
cal fields <5yO Sj a, #s,a (in the co-ordinate representation), 
over which the functional integrals are taken. For the 
sake of simplicity, we write the Hamiltonian (|lll) in the 
phononic limit, where the fluctuation wavelengths are 
long compared to the healing length of the quasicon- 
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\z—z'\/lj 



(9) 



(0 [6]) arc reduced to uj£ 

(expliSjiz') - i6j(z)]) = exp <j - -{[0 B (z') - s (z)} 2 ) - and ^l( k ) ~ c2k ' 2 + AJgn m /Ti, where c = y/gn m /m 

[2 is the speed of sound. Of course, the phase-density de- 

scription can be extended into short-wavelength excita- 
tion range fiol [ll| , bringing about the Hamltonian terms 
oc (dSp s . a /dz) 2 and thus revealing the full Bogoliubov- 
like spectra (O [5]). However, we are not interested in the 
short- wavelength limit, since the respective length scales 
cannot be resolved by optical imaging systems [i, \T3l [I3| . 
The system's description by Eq. (fTTj) is fully consistent 
with Haldane's bosonization method [l6| . The relative 
phase 9 a is accessible through interference patterns ob- 
served in time-of-flight experiments [1, [H, Q • We de- 
velop here the way to evaluate its correlation properties. 
Since the density fluctuations are small, we can decouple 
symmetric and antisymmetric modes [17| and integrate 
out the variables of the symmetric mode. The absence 
of cross-terms containing both Sp a and 9 a in Eq. (fTT|) 
allows us to integrate out Sp a as well and to obtain, as an 
intermediate result, the partition function in the form 



decreases ex exp[— \z — z'|/(2At)] at \z — z'\ — > oo be- 
cause of the unlimited growth of the fluctuations of the 
symmetric component of the phase along the z-direction. 
The correlation properties of the symmetric mode can be 
experimentally measured using the density-density corre- 
lations of the ultracold gas in a time-of-flight experiment 
[l5j . however, this subject is beyond the scope of our 
present paper. 

The phase locking of the relative phase becomes most 
apparent if we treat the evolution of the relative phase 
along z in the harmonic approximation as the Ornstein- 
Uhlenbeck stochastic process Q: while thermal excita- 
tions result in the relative phase diffusion, with the diffu- 
sion coefficient proportional to , the tunnel coupling 
gives rise to the "friction" force that tends to restore a 
small (ultimately zero) local phase difference between the 
two quasicondensates. 



III. CORRELATION FUNCTIONS AND THE 
INTERWELL COHERENCE BEYOND THE 
HARMONIC APPROXIMATION 

A. Equilibrium theory 

In the present work we make a step further with respect 
to the theory of Ref. [l[ and abandon the assumption 
of small phase fluctuations (but still consider small den- 
sity fluctuations, which is a reasonable approximation for 
quasicondensates with repulsive interactions) . We evalu- 



Z = const / T>9 a exp < — dz 



h 2 niu / 80 a 
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that was considered long ago [3 [l9| in the context of the 
statistical mechanics of systems describable by the sine- 
Gordon equation, which is known to adequately account 
for the low-ener gy p hysics of tunnel-coupled ID ultracold 
atomic systems [17l |. 

Note that anharmonic Hamiltonian terms, which de- 
pend on the density fluctuations neglected in our present 
theory, do not affect much the static properties of the 
quasicondensate [ll|. One needs to take them into ac- 
count in the analysis [20| of a slow process of the sys- 
tem's relaxation towards equilibrium starting from a non- 

character- 
_ < T+ for 



equilibrium, pre-thcrmalized initial state [21 
ized by two different temperatures T + and T. 
the symmetric and antisymmetric modes, respectively. 



4 



The applicability range of our fully classical approach 
can be determined as follows. Consider, for the sake of 
simplicity, distances shorter than lj. The effects of tun- 
nel coupling can be neglected at such short length scales, 
and the fully classical correlation function can be esti- 
mated [l| as (exp[i6 a (z')—i6 a (z)]) w cxp(— 2|z — z'\/Xt)- 
We have to compare this result to the power-law de- 
cay of correlations due to quantum effects, which is 
obtained in the limit T -> [HI [ll|. Neglecting, 
as previously, the contribution of the density fluctu- 
ations, we can write lini7'_ ) .o(ex:p[i# a (V) — i6 a {z)\) « 
lim r ^ o ('0l('2')^i(- 2; ))(^l( z )'02(2 / )) and, finally, 



lim (exp[i# a (z') - i0 a (z)]) 

T— >0 



A 



uv 



\z-z'\ 



1/JC 



(13) 



where the quantum mechanical average over the ground 
state is taken, K. = irh^/niT>/(mg) is the Luttinger liq- 
uid parameter (for quasicondensates, which are weakly 
interacting systems, K > 1), and Auv is the ultraviolet 
cutoff of the theory. Eq. (fl3|) is valid if 



\z - z'\ » A 



uv- 



(14) 



The estimation by Popov [22j yields Auv ~ £• 

We can fully neglect quantum fluctuations if their con- 
tribution to the decay of correlations is small, compared 
to the contribution of the thermal noise, on a given length 
scale. The correlation decay is dominated by the thermal 
noise if the classical formula exp(— 2\z — z'\/\t) yields 
stronger decay of correlations than the quantum limit 
(HI, i.e., if 



2\z - z'\/\ T Z /C-Mnflz - ar'l/0 . 



(15) 



The experimentally relevant range of \z — z'\ is bound 
from below by Az opt , as we discussed in Sec. [TT1 and 
Az opt ^> £ in a typical experiment Q . Therefore the use 
of the fully classical approach is reasonable for 



knT> 



mc 



2 £ln(Az opt /£) 



7rAz, 



(16) 



opt 



We can evaluate the partition function (fT^j) using the 
transfer operator technique 0, [ljl |2j|. First of all, we 
evaluate the phase-correlation function as (see Appendix 
lAl for the sketch of derivation) 

(exp[i0 a (z') - i6 a (z)}) = 



OO 

E 

Tt = 



(n|e i9 |0)| \ xp [-(e n -e )\z-z'\} , (17) 



where 



{ne 



'|0) 



(18) 



1 $ n (9) is the eigenfunction (normalized to 1) of the aux- 
iliary Schrddinger-type equation 



At 



(cost 



1 



9 n {0) = e n V n {0), (19) 



and e n , n = 0,1,2,..., is the respective eigenvalue. 
For simplicity, we set periodic (and not quasiperiodic) 
boundary conditions to Eq. ([T9| with the period 27T, 
thus neglecting the band structure of its spectrum, since 
the zero-quasimomentum solutions define all the system 
properties fl9j . which are relevant to our present work. 

In the limit of strong tunnel coupling, lj <C At, the 
operator in the left-hand-side of Eq. (|19[) can be approx- 
imated by the harmonic oscillator Hamiltonian (in proper 
units), and e„ = lj 1 (n + |), n = 0, 1,2, . . . . In this limit 
Eq. (fTT)) reproduces the result © that holds for small 
phase fluctuations. 

In the opposite limit, Eq. (fT9|) can be solved perturba- 
tively, and we obtain 



(cxp[i0 a (z')-i0a(z)]>« -7J + 



Ay 



8iJ 



Xrp 



exp 



2|z-z'| 



l j » A T . (20) 



In what follows, we will be interested in calculating the 
value of 



(cos 9 a 



|cos0|O>, 



(21) 



which can be viewed as the mean interwell coherence. 
This expression can be derived in different ways, e.g., 
from Eq. (|17|) by employing the statistical independence 
of phase fluctuations at two very distant points, \z— z'\ — > 
oo, and recalling that (sin# a ) =0. In a general case, Eq. 
(|2T|) can be evaluated from the lowest-energy solution of 
the Mathieu equation [24| . In the two limiting cases we 
obtain the asymptotics 



exp(-Zj/Ar), h < A T 
" /(8l 2 j), lj » At 



A j ■ 



(22) 



A possible physical explanation of the fact that the mean 
interwell coherence decreases at IjJXt — > oo much slower 
that the harmonic approximation [1( predicts, is the large 
probability of thermal excitation of a soliton in this limit. 
Each emerging soliton decreases the number of phononic 
states by 1 [T{|, and the phononic density of states is 
reduced mostly in the long- wavelength range (for phoiion 
momenta less than or of the order of U/lj), which gives 
the main contribution to the long-distance behavior of 
the correlation function ([T71) and, hence, to (cos# a ). 



B. Relaxation to the equilibrium after a quench 



The results of Sec. MI Al are obtained at the equilib- 
rium. However, it is interesting to investigate also the 
process of equilibration in the system of two ID quasi- 
condensates after a quench. The study of this dynami- 
cal problem is motivated by our recent numerical results 
[25j related to thermalization in a single ID quasicon- 
densate. In Ref. (2f| we found that, despite the nu- 
merically confirmed intcgrability of the system, phononic 
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(low-momentum) modes rapidly relaxed from their ini- 
tial non-equilibrium state towards a final equilibrium 
state; particle-like (large-momentum) excitations, on the 
contrary, exhibited almost no relaxation. The equilib- 
rium ensemble of phonons was different from the classi- 
cal limit of equipartition of the thermal energy between 
all the degrees of freedom and was quite close to the 
Bosc-Einstcin distribution with the temperature T e g de- 
termined by the total excitation energy of the initial 
non-equilibrium state. Observed fluctuations around this 
equilibrium state were due to the finite size of the system 
inherent to numerical modeling. Remarkably, the cor- 
relations observed at the length scales, which are large 
compared to the healing length to the healing length, 
as well as to the wavelength of an elementary excitation 
with the energy equal to /ceT^, were well described by 
classical expressions. Note that the main contribution 
to the noise on these length scales stems from the low- 
energy excitations, which approximately exhibit classical 
equipartition of energy. 

The need to extend the numerical approach of Rcf. 
[2|| to tunnel-coupled ID quasicondensates can also be 
seen from the following considerations. Our aim is to nu- 
merically check the theoretically predicted correlations 
of two tunnel-coupled quasicondensates at equilibrium. 
This equilibrium state can be viewed as a result of the 
system's relaxation from its initial non-equilibrium state. 
Moreover, the available analytic theory predicts only av- 
erages; unlike the case of harmonic approximation, there 
is no way yet to generate individual realizations of the 
phase, obeying the necessary statistics, without simulat- 
ing numerically the equilibration process. The most ob- 
vious way to obtain numerically the equilibrium solution 
is to observe the numerical relaxation after a quench and 
wait until a steady-state regime establishes. Particular 
type of the quench and the corresponding initial condi- 
tions are, up to a certain degree, arbitrary, as long as the 
system exhibits true relaxational dynamics. 

Motivated by these considerations, we performed nu- 
merical modeling of the thermal equilibrium values of 
(cos 9 a ) after the dynamical process of relaxation in our 
system after a quench. We simulated the time evolu- 
tion of two coupled Gross-Pitaevskii equations using the 
split-step method (2(| previously used by us [25[ to sim- 
ulate the dynamics of a single quasicondensate and now 
extended to the case of tunnel-coupled systems. As the 
initial conditions we took two independent quasiconden- 
sates with phonon modes populated randomly according 
to the Bose-Einstein thermal distribution. At t = we 
quenched the system by switching on the tunnel coupling 
between them. We solved this coupled system for a time 
long enough to provide equilibration. 

To juxtapose the input parameters of our numerical 
simulations to typical parameters of modern atom-chip 
experiments [1, OS [l4[ , we give the system parameters 
used in our simulations first in dimensional units, but 
later show them also in dimensionlcss form. The lin- 
ear density for a single quasicondensate tt-id = 30 jimT 1 




FIG. 1: (Color online) Mean interwell contrast as a function of 
the ratio of the length scales At and lj. Solid line: exact the- 
ory given by Eq. ([21]). Dashed line: small- fluctuations approx- 
imation (cos 6a) — exp(— Ij/Xt) following from the linearized 
theory Dots: results of the numerical simulations of the 
equilibration dynamics of two coupled condensates. Units on 
the axes are dimensionless. Inset: Magnified part of the main 
plot for small \t/Ij, illustrating the high-temperature asymp- 
totics of Eq. (|21[) in comparison to the linearized theory result. 



and the interaction constant g = 2 7kjj_a s with the ra- 
dial trapping frequency = 2ir x 3 kHz and the s-wave 
scattering length a s = 5.3 nm for 87 Rb yields the healing 
length £ w 0.35 fxxa and the Luttinger liquid parame- 
ter K. w 33. The periodic boundary conditions were set 
at an interval of the length L = 100 /Ltm ps 290 £. The 
maximum integration time was t max = 0.8 s. After few 
hundreds milliseconds some kind of equilibrium was ob- 
tained. The total energy of the system was conserved 
in our numerical simulations with a good (~ 10 -3 ) ac- 
curacy, however, it was constantly redistributed in an 
oscillatory manner between different low-frequency ele- 
mentary modes, including Josephson oscillations. The 
nonlinear interaction between different modes (see Sec- 
tion HVj) lead to excitation of Josephson oscillations of 
the total number imbalance (Ni — A^)/2, where Nj is 
the integral of the density in the jth quasicondensate 
over the whole length L, i.e., the number of atoms in 
this quasicondensate, Nx + N^ = 2N. In general, the nu- 
merical stability of our split-step method was controlled 
using the criteria of Ref. [27| . The thermal coherence 
length was determined from the phase-correlation func- 
tions for each of the two quasicondensates taken sepa- 
rately by comparison of the numerically obtained value of 
(cxp[i9j(z) — i6j(z')]}, j = 1,2, with its theoretical value 
cxp(— \z — z'\/Xt) for \z — z!\ < lj [HHHl (if we trace out 
the phase and density variables of one of the two tunnel- 
coupled quasicondensates, the properties of its remaining 
counterpart will be described by the same temperature 
as of the whole system at equilibrium). The averaging is 
performed over statistically uncorrelatcd (separated by 
sufficiently large distances) intervals of the whole length 
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L for \z 



< At- We never obtain complete equili- 



bration. In each realization, the correlation length At 
obtained in such a way oscillates around certain mean 
value, and so does the value of (cos# a ) (averaged over 
the length L). Typically, At ~ 8 /im, which corresponds 
to T ps 40 nK. 

We present the results of our numerical simulations in 
Fig-ffl Dots represent mean values of (cos 9 a ) obtained by 
averaging over both the time (on the quasi-cquilibration 
stage of the system evolution) and the ensemble of real- 
izations. The error bars in Fig.[T]show the standard devi- 
ations of (cos# a ) and At- These error bars indicate slow, 
quasipcriodic variations of (cos 9 a ) and At detected in 
our simulations. The range of Xt/Ij shown in Fig. [1] cor- 
responds to J increasing from 2ir x 0.1 Hz up to 2tt x 8 Hz. 

To summarize the results of the present Section, we 
can state that we developed a theory describing the static 
correlation properties more precisely than the harmonic 
model [l| . Our approach is based on consideration of the 
classical partition function for the antisymmetric mode 
of our problem (describablc by the sine-Gordon model) 
and application of the well-known transfer operator tech- 
nique [ill, [l9| . As one can see from Fig. [TJ the difference 
between our results and those of Ref. [l| is most apparent 
for intermediate and small values of At /lj (intermediate 
and weak tunnel coupling). 



IV. JOSEPHSON OSCILLATIONS IN A NOISY 
EXTENDED JUNCTION 
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FIG. 2: (Color online) Josephson oscillations for J = 27rx8 Hz 
(for other system parameters see Section IIIII of the main 
text), (a) The relative imbalance as a function of time, (b) 
The power spectrum of the atom-number imbalance (aver- 
aged over 7 realizations), peaked at theoretically predicted 
ljj/(2tz) = 157 Hz and broadened by thermal fluctuations. 



The thermal noise effects considered in Sec. IIIII reduce 
the frequency of Josephson oscillations. 

Consider the absolute number imbalance between two 
wells, N r2 = (Ni — AT 2 )/2, and its canonically conjugate 
variable, the overall phase difference 4> between two qua- 
sicondensates. In the limit of the atomic repulsion energy 
dominating over the tunneling, gniD = gN/L 3> hJ, and 
for small-amplitude oscillations, \N± — N%\ <C N, the evo- 
lution of these "global" variables is described by the set 
of equations (see Appendix |B|) 



2gN r 
Lh 



d 
dt 



—N12 = 2Jn 1D / dz s'md a , 
at J 



(23) 



(24) 



which is reduced, after elimination of the number- 
difference variable, to 



is the frequency of the Josephson oscillations for bosonic 
junction unaffected by thermal noise. At zero tempera- 
ture, when the thermal noise is absent, and for ln(L/£) <C 
/C, when the quantum noise can be neglected, spatial ex- 
tension of the ultracold-atomic Josephson junction plays 
no role and we can derive Eq. ([26]) from the results of Ref. 
Q . In the case of small-amplitude Josephson oscillations, 
the statistical properties of cos# a and cos(# a — $) do not 
differ significantly, in particular, (cos 6* a ) ps (cos(6* a — $)), 
i.e., the quadratic in $ correction is negligible, and Eq. 
^5l) reduces to 



where 



— $ + [ w 3+*w3(t)]$ = c(t), 



w .i = wj (cos6> a ). 



(27) 



(28) 



In Eq. (|27[) we explicitly indicate the time argument 
of the random driving force 



where 



d 2 1 f L 

^ = -^°lJ dzM 



(25) 



(26) 



and the term 



'LJ 



C(t) =wj - / dz sm(9 a -$) 



(29) 



6wj(t) = wj — / dz (cos6» a - (cos0 a )) (30) 
L Jo 
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FIG. 3: (Color online) The same as in Fig. [2] but for J = 2n x 
0.1 Hz (irregular behavior). The spectral peak at theoretically 
predicted ujj/(2n) = 6.5 Hz is smeared out. S 7^ at uj = 
due to finite integration time. 
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FIG. 4: (Color online) Dots: the square of the frequency of 
numerically obtained Josephson oscillations (normalized to its 
value oijq for zero thermal noise) as a function of the mean 



interwell coherence (cos# a )- 
prediction given by Eq. ([28 
sionless. 



The straight line: theoretical 
Units on the axes are dimen- 



that describes fluctuations of the oscillation frequency 
due to the noise of 9 a caused by excitations with non- 
zero momenta. 

If we prepare some appreciable initial imbalance at t = 
0, we obtain, to the first approximation, free Josephson 
oscillations governed by the equation d 2 $>/dt 2 + ui 2 $ = 0, 



i.e., with the frequency reduced by \J (cos a ) compared 
to the noise- free case of Eq. ((25)1. 

The presence of the noise broadens the power spectrum 
of Josephson oscillations 



S(u) 



dte lut r](t) 



(31) 



relative number 
(|3ip is taken 
system has al- 
state (typically, 
the theory |28j 



where 77 = (N x - N 2 )/{2N) is the 
imbalance. The integration in Eq. 
over the time interval r when the 
ready reached its nearly-equilibrium 
r w 0.65 s). If is high enough, 
predicts S(lu) to be a peaked function, centered at 
wj and having the half-width at the half-maximum of 
the peak height 7 = [hL/{8gk B T)}Rc dt'(((t)((t + 
i')) exp(iwjt'). The latter expression, roughly evaluated 
as 7 ~ ^k&T/ (%K.(cos a ) 2 ), correctly describes the order 
of magnitude of the bandwidth Alo/(2tt) ~ 10 Hz of the 
numerically obtained spectra S(uj). 

The presence of the random driving force is the source 
of excitation of Josephson oscillations in the course of 
the system's evolution, even if initially at t = 0, $ = 
and 77 oc ^<f> = 0. Note, that all the elementary excita- 
tions with nonzero momenta in the antisymmetric mode 
have frequencies larger than than ojj. The energy trans- 
fer between nonzero-momentum excitations and Joseph- 
son mode is thus an essentially nonlinear process. The 
nonlinear structure of the right-nand-side of Eq. ((29)) 
provides the presence of the frequency wj in the spec- 
trum dt'(((t)C{t + t')) exp(iwt') of the driving force 
and thus ensures the parametric excitation of the Joseph- 
son oscillations. 

We confirmed our analytic estimations by the numeri- 
cal simulations of two coupled ID Gross-Pitacvskii equa- 
tions already described in Subsection IIIIBI An example 
of a sharp-peaked power spectrum of relative number im- 
balance is given in Fig. [2J together with an example of 
time dependence of r]. 

If, on the contrary, wj <C wt, where lot = 2c/Ar is 
the typical time scale of fluctuations of £(<), then the 
behavior of rj(t) becomes irregular and S(u>) does not 
exhibit a peak at u> ~ wj any more (see Fig. [3]). 

The results of numerical simulations shown in Figs. 
[2] and [3] demonstrate certain energy exchange, but no 
full equilibration between the Josephson oscillations and 
phononic modes. If we set 3>|t=o = and rj\ t= o = 
for J/(2tt) = 8 Hz (or 0.1 Hz), then at times t between 
650 ms and 1 s the mean energy of Josephson oscillations 
is by an order of magnitude (or by 1.5 orders of magni- 
tude, respectively) less than ks,T, where temperature T 
is determined from the phase-correlation function for a 
single quasicondensate and is thus associated with the 
phononic modes. This may indicate an extremely long 
thcrmalization time for Josephson oscillations. 

We selected our simulations that display a pronounced 
narrow peak of S(uj) far from zero frequency (which was 
the case for J > 2tt x 0.7 Hz), estimated the Josephson 
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frequency uij and analyzed the dependence of Wj on the 
mean interwell coherence. The resulting values are in a 
good agreement with our theoretical prediction given by 
Eq. ([25)1 , as can be seen from Fig. 01 



V. CONCLUSION 

To conclude, we applied the transfer-operator tech- 
nique to evaluate coherence and correlation properties of 
two tunnel-coupled ID weakly-interacting, ultracold sys- 
tems (quasicondensates) of bosonic atoms. These proper- 
ties are determined by the ratio of the two length scales: 
At that describes the spatial scale of the loss of correla- 
tions between two points and lj that describes the scale 
for the phase-locking between two quasicondensates due 
to interwell tunneling. In the limit lj < At the fluctua- 
tions of the relative phase are small and we reproduce the 
results of the linearized theory of Ref. [lj . In the opposite 
case, we found the mean interwell coherence to decrease 
much slower (oc A T /Z^) than the exponential law pre- 
dicted by the linearized theory. We interprete such a be- 
havior as a signature of thermal creation of sine-Gordon 
solitons, which provide a shift of the relative phase by 
2w and thus do not contribute to the coherence loss, and 
the corresponding decrease of the density of states for 
phonons (the excitations responsible for the coherence 
loss at large distances). 

Our analytic estimations are confirmed by numerical 
modeling of the equilibrium state as a final state of the 
system's relaxational evolution after a quench. This task 
is solved by extending our numerical method [25| to in- 
tegration of two coupled ID Gross-Pitaevskii equations. 

We demonstrate, both analytically and numerically, 
that thermal fluctuations of the relative phase between 
two quasicondensates reduce the frequency of Joscphson 
oscillations in proportion to y/ (cos 9 a ) and broaden their 
spectrum. If the theoretically predicted value of loj is 
much less than the bandwidth of the thermal fluctuation 
(which is of the order of the speed of sound divided by 
At), regular Josephson oscillations are not observed. 

This work was supported by the the FWF (Project 
No. P22590-N16). The authors thank T. Berrada and J. 
Schmiedmayer for helpful discussions. 



Appendix A: Derivation of Eq. (|17[) 

We briefly recall here the basics of the transfer operator 
technique, following Refs. [H, [l^,[2^. We introduce a 
lattice with the step Az = L/M, M being the number of 
sites. We assume cyclic boundary conditions, 

QaM+l = #al- (Al) 

Then the partition function (|12[) can be written as 
Z = J d9 a i . . . J dd a M J d9 a M+l <5(#aM+l — #al) x 



]Jexp[-f(e aj ,e aj+1 )}, 

3=1 



(A2) 



where 

f(O a j, #aj+l) 



Amk B TAz v a] 
TiJniYyAz 



(2 - COS0 a j - COS0 a j + i) (A3) 



and integrals in our case are taken from — it to ir. We 
omit the constant prcfactor in Eq. (|A2j) for the sake 
of simplicity. Assume that eigenfunctions ^ n (9) of the 
transfer operator 

d9 aj e-W**< e *^* n (6 aj ) = e-^ Az ^ n (9 aj+1 ) 

(A4) 

form a set, which is complete, orthogonal, and normalized 
to unity, namely 



d0*;,(0)* B (0) 



5>;(*o*»(0) = 5(9' 



(A5) 
(A6) 



Substituting Eq. (|A6|) into Eq. (|A2|) and using Eq. (|A4|) . 
we obtain 



Z = y^cxp(-e„L). 



(A7) 



The eigenvalues e n are positive; in the thermodynamic 
limit the partition function (|A7|) is dominated by the 
lowest eigenvalue tQ, 



Z k, exp(-e £), 



L 



oo. 



(A8) 



In the continuous limit Az — > Eq. (|A4[) is equivalent 
to the Schrodinger-type equation (|T9| . Strictly speaking, 
the spectrum of Eq. (|A4[) is shifted with respect to the 
spectrum of Eq. (|19p by a common offset so, which is 
related to normalization of the eigenfunctions. Since Sq 
does not depend on n, we neglect it in our calculations. 

To calculate correlation functions, in particular, Eq. 
CI]), we note that e l6 ^ z ") and e~ l6 ^ z) act on * like 
quantum-mechanical perturbations, coupling ^/q to the 
whole spectrum of eigenfunctions with the matrix ele- 
ments given by Eq. (|18[) . Therefore the leading term for 
(cxp[i9 a (z')—i9 a (z)]) in the limit of L — > oo is the second- 
order pcrturbative correction to the propagator for the 
ground state (with L playing the role of imaginary time), 
and we obtain thus Eq. p7)l . 



Appendix B: Derivation of Eqs. (f23l |24|) 

We begin with the lattice version of the classical sine- 
Gordon Hamiltonian that describes the dynamics of the 



antisymmetric mode of our system: 

,2, 



M 

^ = E 

3 = 1 



h niD 



2hJnn)Az(l — cos# a:; -) 



(Bl) 



,Az is 



where the jth generalized co-ordinate SN a j - 
the half-difference of the atomic numbers in the 1st and 
2nd quasicondensates at the jth site, i.e., the variable 
canonically conjugate to the local phase difference 9 a j 
(the jth generalized momentum). Here we neglect the 
nonlinear coupling between the symmetric and antisym- 
metric modes, like in Eq. (|lll) in the continuous limit. 

For the sake of simplicity, we assume an odd number of 
sites in the lattice, M = 2Mq + 1, where Mq is a positive 
integer. Then we do a canonical transformation 



Mo 

w aj = £ SN^ewtj), 9 aj = J2 WvVJ), 



t=-M 



M 

E 

C=-Ma 



where 



(B2) 



V^/Mcos(27r£j/M), I = -1, -2, . . . , -M 
= { 1/VM, £ = 



(B3) 



y/2/M sin(2Tr£ j/M), £ = 1,2,...,M 



Then the Hamiltonian (IB1I) reads 



H a = 




[l-cos(2^/M)]0lW 



(B4) 



M 



2h Jriiv Az < 1 — cos 

3 = 1 I 

From the Hamiltonian equations 



Mo 

E Wv(t,j) 

£=-M 



at d9 a (£) at 



we find, in particular, 



dH a 



dSN a (£) 



(B5) 



(B6) 



d 
rft 



A I 



sin 

j'=i 



' M 

E Ut)n{z,3) 

i=-M 



(B7) 

In the limit of Az — > the sums over j converge to inte- 
grals over z. Taking into account that N12 = J dz5p a = 



^2j = i5N a j = \ / AlSN a (0), identifying the generalized 
momentum conjugate to N\ 2 as <I> = 9 a (0)/y/M = 
{l/M)Y J f =x O a j and recalling that L = MAz, we ob- 
tain Eqs. ([231 CM]) . The spatially fluctuating part of the 
phase is then 9 a — $. 



[1] N. K. Whitlock and I. Bouchoule, Phys. Rev. A 68, 

053609 (2003). 
[2] I. Bouchoule, Eur. Phys. J. D 35, 147 (2005). 
[3] I. Lesanovsky and W. von Klitzing, Phys. Rev. Lett. 98, 

050401 (2007). 

[4] R. Hipolito and A. Polkovnikov, Phys. Rev. A 81, 013621 
(2010). 

[5] H.-P. Stimming, N. J. Mauser, J. Schmiedmayer, and I. 

E. Mazets, Phys. Rev. Lett. 105, 015301 (2010). 
[6] J. Brand, T. J. Haigh, and U. Ziilicke, Phys. Rev. A 81, 

025602 (2010). 

[7] T. W. A. Montgomery, R. G. Scott, I. Lesanovsky, and 

T. M. Fromhold, Phys. Rev. A 81, 063611 (2010). 
[8] T. Betz, S. Manz, R. Biicker, T. Berrada, Ch. Roller, 

G. Kazakov, I. E. Mazets, H.-P. Stimming, A. Perrin, T. 

Schumm, and J. Schmiedmayer, Phys. Rev. Lett. 106, 

020407 (2011). 

[9] A. Smerzi, S. Fantoni, S. Giovanazzi, and S. R. Shenoy, 
Phys. Rev. Lett. 79, 4950 (1997). 
[10] V. N. Popov, Functional Integrals and Collective Excita- 
tions (Cambridge University Press, Cambridge, 1987). 
[11] C. Mora and Y. Castin, Phys. Rev. A 67, 053615 (2003). 
[12] M. A. Cazalilla, J. Phys. B 37, SI (2004). 



[13] S. Hofferberth, I. Lesanovsky, B. Fischer, T. Schumm, 
and J. Schmiedmayer, Nature 449, 324 (2007). 

[14] S. Hofferberth, I. Lesanovsky, T. Schumm, A. Imam- 
bekov, V. Gritsev, E. Demler, and J. Schmiedmayer, Na- 
ture Phys. 4, 489 (2008). 

[15] S. Manz, R. Biicker, T. Betz, Ch. Roller, S. Hoffer- 
berth, I. E. Mazets, A. Imambekov, E. Demler, A. Per- 
rin, J. Schmiedmayer, and T. Schumm, Phys. Rev. A 81, 
031610(R) (2010). 

[16] F. D. M. Haldane, Phys. Rev. Lett. 47, 1840 (1981). 

[17] V. Gritsev, A. Polkovnikov, and E. Demler, Phys. Rev. 
B 75, 174511 (2007). 

[18] J. A. Krumhansl and J. R. Schrieffer, Phys. Rev. B 11, 
3535 (1975). 

[19] J. F. Currie, J. A. Krumhansl, A. R. Bishop, and S. E. 

Trullinger, Phys. Rev. B 22, 477 (1980). 
[20] H.-P. Stimming, N. J. Mauser, J. Schmiedmayer, and I. 

E. Mazets, Phys. Rev. A 83, 023618 (2011). 
[21] M. Gring, M. Kuhnert, T. Langen, T. Kitagawa, B. 

Rauer, M. Schreitl, I. Mazets, D. Adu Smith, E. Demler, 

and J. Schmiedmayer, Science 337, 1318 (2012). 
[22] V. N. Popov, Pis'ma Zh. Eksp. Teor. Fiz. 31, 560 (1980) 

[Sov. Phys. JETP Lett. 31, 526 (1980)]; see also [lj. 



10 



[23] D. J. Scalapino, M. Sears, and R. S. Ferrell, Phys. Rev. 

B 6, 3409 (1972). 
[24] M. Abramowitz and I. A. Stegun, eds., Handbook of 

Mathematical Functions (Dover, NY, 1972), Ch. 20. 
[25] P. Grisins and I. E. Mazets, Phys. Rev. A 84, 053635 

(2011). 

[26] M. Thalhammer, M. Caliari and C. Neuhauser, J. Com- 



put. Phys. 228, 3 (2009). 
[27] J. A. C. Weideman and B. M. Herbst, SIAM J. Numer. 

Anal. 23, 485 (1986). 
[28] M. Gitterman, The Noisy Oscillator (World Scientific, 

Singapore, 2005), Ch. 8. 



